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^ ■ 

o : 

' Abstract. We define quasiconvex programming, a form of generalized linear programming in which 
one seeks the point minimizing the pointwise maximum of a collection of quasiconvex functions. We 

' survey algorithms for solving quasiconvex programs either numerically or via generalizations of the 

1^ , dual simplex method from linear programming, and describe varied applications of this geometric opti- 

■ mization technique in meshing, scientific computation, information visualization, automated algorithm 
' analysis, and robust statistics. 

I— i! 1 Introduction 

O: 

■ Quasiconvex programming is a form of geometric optimization, introduced by Amenta et al. in 
c/3 , the context of mesh improvement techniques [3] and since apphed to other problems in meshing, 

I ^/ scientific computation, information visuahzation, automated algorithm analysis, and robust statis- 
tics [8,9, 14,33]. If a problem can be formulated as a quasiconvex program of bounded dimension, it 

^ . can be solved algorithmically in a linear number of constant-complexity primitive operations by gen- 

\Q I eralized linear programming techniques, or numerically by generalized gradient descent techniques. 

' In this paper we survey quasiconvex programming algorithms and applications. 

<N : 

' 1.1 Quasiconvex Functions 

o ■ 

' Let y be a totally ordered set, for instance the real numbers M or integers Z ordered numerically. 

Q . For any function f : X Y, and any value X €Y, we define the lower level set 

-t.: f^^ = {xex\f{x)<x}. 

rS ; 

■ A function q : X i-^ Y, where X is a convex subset of M°', is called quasiconvex [23] when its 
lower level sets are all convex. A one-dimensional quasiconvex function is more commonly called 
unimodal, and another way to define a quasiconvex function is that it is unimodal along any line 
through its domain. 

As an example, let H = {{x,y) [ y > 0} be the upper halfplane in M^, let u = (—1,0) and 
w = (1,0), and let q measure the angle complementary to the one subtended by segment uw from 
point v: q{v) = 180° — Zuvw. Then, each level set q-^ consists of the intersection with H of a disk 
having u and w on its boundary (Figure^. Since these sets are all convex, q is quasiconvex. 

Quasiconvex functions are a generalization of the well-known set of convex functions, which are 
the functions M'^ M satisfying the inequality f{px + (1 — p)y) < p f{x) + (1 — p)f{y) for all 
x,y € M'^ and all < p < 1: it is a simple consequence of this inequality that any convex function 
has convex lower level sets. However, there are many functions that are quasiconvex but not convex; 
for instance, the complementary angle function q defined above is not convex, as can be seen from 
the fact that its values are upper bounded by 180. As another example, the function xk{x) that 



(-1,0) 



(1,0) 



Fig. 1. Level sets of the quasiconvex function q(v) = 180° — Zuvw, for u = (—1,0) and w = (1,0), 
restricted to the halfplane y > 0. 

takes the value within a convex set K and 1 outside K has as its lower level sets K and R'^, so 
is quasiconvex, but not convex. 

If r is convex or quasiconvex and f :¥ Z is monotonically nondecreasing, then q{x) = f{r{x)) 
is quasiconvex; for instance the function xk above can be factored in this way into the composition 
of a convex function dxix) measuring the Euclidean distance from xto K with a monotonic function 
/ mapping to itself and all larger values to 1. In the other direction, given a quasiconvex function 
q : X ^ Y , one can often find a monotonic function / : 1" i— > M that, when composed with g, 
turns it into a convex function. However this sort of convex composition is not always possible. For 
instance, in the case of the step function xk described above, any nonconstant composition of xk 
remains two-valued and hence cannot be convex. 

1.2 Nested Convex Families 

Quasiconvex functions are closely related to nested convex families. Following Amenta et al. [3] , 

we define a nested convex family to be a map k : Y ^ K^W^), where y is a totally ordered set 
and K{W^) denotes the family of compact convex subsets of W^, and where n is further required to 
satisfy the following two axiomatic requirements (the second of which is a slight generalization of 
the original definition of Amenta et al., that allows Y to be discrete): 

1. For every Ai,A2 G Y with Ai < A2, k(Ai) C k{X2)- 

2. For all A G y for which A = inf{A' | A' > A}, k{\) = {^^^^ k{\'). 

If Y has the property that every subset of Y has an infimum (for instance, Y 
then from any nested convex family k :Y ^ K{W^) we can define a function q,^ 
following formula: 

(Ik{x) = inf { A I X G k{X) } . 

Lemma 1. For any nested convex family k :Y t-^ K{W^) and any X eY, q^^ = k{X). 

Proof. The lower level sets of are 

q^^ = {xeR'^\ q^ix) < X}= {xeR'^linf {X' \xe k{X') } < A}. 

For any x G k(X), A G { A' | x G k{X') } so the infimum of this set can not be greater than A and 
X G q^^. For any x ^ k(A), inf { A' | x G k(A') } > A"*" > A by the second property of nested convex 
families, so x ^ q^^. Therefore, q^^ = k{X). □ 
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In particular, has convex lower level sets and so is quasiconvex. 

Conversely, suppose that q is quasiconvex and has bounded lower level sets. Then we can define 
a nested convex family 



f]y^^cl{q^^') if A = inf{A' I A' > A} 
cl{q-^) otherwise 



where cl denotes the topological closure operation. 

If q does not have bounded lower level sets, we can still form a nested convex family by restricting 
our attention to a compact convex subdomain K C M'^: 



Kg,ii-(A) 



nA'>A cl(i^ n q^^') if A = inf{A' | A' > A} 
cl{K n q-'^) otherwise 



This restriction to a compact subdomain is necessary to handle linear functions and other functions 
without bounded level sets within our mathematical framework. 

The following two theorems allow us to use nested convex families and quasiconvex functions 
interchangeably for each other for most purposes: more specifically, a nested convex family conveys 
exactly the same information as a continuous quasiconvex function with bounded lower level sets. 
Thus, later, we will use whichever of the two notions is more convenient for the purposes at hand, 
using these theorems to replace an object of one type for an object of the other in any algorithms 
or lemmas needed for our results. 

Theorem 1. For any nested convex family 

Proof. If A is not an infimum of larger values, then q^ix) < A if and only if x G k(A). So Kq„(A) = 
cl(g^^^) = {x I q^ix) <X} = k{X). 
Otherwise, by Lemma ^ 

«..(A)= n ^i('^(A')) 

A'>A 

The closure operation does not modify the set k(A'), because it is already closed, so we can replace 
c1(k(A')) above by k{X')), giving 

A'>A 

The intersection on the right hand side of the equation further simplifies to k{X) by the second 
property of nested convex families. □ 

Theorem 2. If q : X M. is a continuous quasiconvex function with bounded lower level sets, then 
qK, = q- 

Proof. By Lemma^ q^^ = Kq{X). Assume first that A = inf{A' | A' > A}. Expanding the definition 
of Kg, we get 

A'>A 

If q is continuous, its level sets are closed, so we can simplify this to 



A'>A 
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Suppose the intersection on the right hand side of the formula is nonempty, and let x be any point 
in this intersection. We wish to show that q{x) < A, so suppose for a contradiction that q{x) > A. 
But then there is a value A' strictly between A and q{x) (else A would not be the infimum of all 
greater values), and x ^ q-^ , contradicting the assumption that x is in the intersection. Therefore, 
q{x) must be at most equal to A. 

As we have now shown that q{x) < A for any x in q^^, it follows that q^^ can not contain any 
points outside q-^. On the other hand, q^^ is formed by intersecting a collection of supersets of 
q-''^, so it contains all points inside q-^. Therefore, the two sets are equal. 

If A 7^ inf{A' > A}, the same equality can be seen even more simply to be true, since we have 
no intersection operation to eliminate. Since g^, ^ind q have the same level sets, they are the same 
function. □ 

Due to these two theorems, we do not lose any information by using the function q,^ in place 
of the nested convex family k, or by using the nested convex family Kq^ = n in place of a quasi- 
convex function that is of the form q = q,^ or in place of a continous quasiconvex function with 
bounded lower level sets. In most situations quasiconvex functions and nested convex families can 
be treated as equivalent and interchangeable: if we are given a quasiconvex function q and need 
a nested convex family, we can use the family Kg, and if we are given a nested convex family k 
and need a quasiconvex function, we can use the function or qK,K- Our quasiconvex programs' 
formal definition will involve inputs that are nested convex families only, but in our applications of 
quasiconvex programming we will describe inputs that are quasiconvex functions, and which will 
be assumed to be converted to nested convex families as described above. 



1.3 Quasiconvex Programs 

If a finite set of functions qi are all quasiconvex and have the same domain and range, then the 
function Q{x) = maxjg5 qi{x) is also quasiconvex, and it becomes of interest to find a point where 
Q achieves its minimum value. For instance, in Section 12.21 below we discuss in more detail the 
smallest enclosing ball problem, which can be defined by a finite set of functions qi, each of which 
measures the distance to an input site; the minimum of Q marks the center of the smallest enclosing 
ball of the sites. Informally, we use quasiconvex programming to describe this search for the point 
minimizing the pointwise maximum of a finite set of quasiconvex functions^. 

More formally. Amenta et al. [3] originally defined a quasiconvex program to be formed by a set 
of nested convex families 5 = K2, ■ ■ ■ the task to be solved is finding the value 



A{S) = inf I (A,x) X G f] 



where the infimum is taken in the lexicographic ordering, first by A and then by the coordinates of x. 
However, we can simplify the infimum operation in this definition by replacing it with a minimum; 
that is, it is always true that the set defined on the right hand side of the definition has a least 
point A(S). To prove this, suppose that {X,x) is the infimum, that is, there is a sequence of pairs 
{Xj,Xj) in the right hand side intersection that converges to (A, x), and {X,x) is the smallest pair 
with this property. Clearly, each Xj > X (else {Xj,Xj) would be a better solution) and it follows 

^ The term quasiconvex programming has also been apphed to the problem of minimizing a single quasiconvex function 
over a convex domain; e.g., see [54,93]. The two formulations are easily converted to each other using the ideas 
described in Section lT^ For the applications described in this survey, we prefer the formulation involving minimizing 
the pointwise maximum of multiple quasiconvex functions, as it places greater emphasis on combinatorial algorithms 
and less on numerical optimization 
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from the fact that the sets Kj are closed and nested that we can take each Xj = x. But then, it 
follows from the second property of nested convex families that x G Ki{\) for all Ki G S. 

In terms of the quasiconvex functions defining a quasiconvex program, we would like to say that 
the value of the program consists of a pair (A,x) such that, for each input function qi, qi{x) < A, 
and that no other pair with the same property has a smaller value of A. However, maxiqi{x) may 
not equal A if at least one of the input quasiconvex functions is discontinuous. For instance, consider 
a one-dimensional quasiconvex program with two functions qo{x) = \x\, qi{x) = 1 for x > 0, and 
qi{x) = for X < 0. This program has value (0, 0), but max{go(0), qi{0)} = 1- The most we can say 
in general is that there exists a sequence of points xj converging to x with limj^oo maxj qi{xj) = A. 
This technicality is, however, not generally a problem in our applications. 

In subsequent sections we explore various examples of quasiconvex programs, algorithms for 
quasiconvex programming, and applications of those algorithms. 

2 Examples 

We begin our study of quasiconvex programming by going through some simple examples of geo- 
metric optimization problems, and showing how they may be formulated as low-dimensional qua- 
siconvex programs. 

2.1 Sighting point 

When we introduced the definition of quasiconvex functions, we used as an example the complemen- 
tary angle subtended by a line segment from a point: q{v) = 180° — Zuvw. If we have a collection of 
line segments forming a star-shaped polygon, and form a quasiconvex program from the functions 
corresponding to each line segment, then the point v that minimizes the maximum function value 
must lie in the kernel of the polygon. If we define the angular resolution of the polygon from v to 
be the minimum angle formed by any two consecutive vertices as seen from v, then this choice of 
V makes the angular resolution be as large as possible. 

This problem of maximizing the angular resolution was used by Matousek et al. [66] as an 
example of an LP-type problem that does not form a convex program. It can also be viewed as a 
special case of the mesh smoothing application described below in Section 14.11 

Earlier, McKay [67] had asked about a similar problem in which one wishes to choose a viewpoint 
maximizing the angular resolution of an unordered set of points that is not connected into a star- 
shaped polygon. However, it does not seem possible to form a quasiconvex program from this version 
of the problem: for star-shaped polygons, we know on which side of each line segment the optimal 
point must lie, so we can use quasiconvex functions with level sets that are intersections of disks 
and halfplanes, but for point sets, without knowing where the viewpoint lies with respect to the 
line through any pair of points, we need to use the absolute value 1^(1^)1 of the angle formed at v by 
each pair of points. This modification leads to non-quasiconvex functions with level sets that are 
unions or intersections of two disks. It remains open whether an efficient algorithm for McKay's 
sighting point problem exists. 

2.2 Smallest Enclosing Ball 

Consider the problem of finding the minimum radius Euclidean sphere that encloses all of a set 
of points S = {pi} C M'^ (Figure 121 left). As we show below, this smallest enclosing ball problem 
can easily be formulated as a quasiconvex program. The smallest enclosing ball problem has been 
well studied and linear time algorithms are known in any fixed dimension [26,35,41,68,92], so the 
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Fig. 2. Smallest enclosing ball of a set of points (left), and the level sets of maxjgj(x) for the 
distance functions qi defining the quasiconvex program for the smallest enclosing ball (right). 

quasiconvex programming formulation does not lead to improved solutions for this problem, but 
it provides an illuminating example of how to find such a formulation more generally, and in later 
sections we will use the smallest enclosing ball example to illustrate our quasiconvex programming 
algorithms. 

Define the function qi{x) = d{x,pi) where d is the Euclidean distance. Then the level set q^^ is 
simply a Euclidean ball of radius A centered at pi, so qi is quasiconvex (in fact, it is convex). The 
function qs{x) = maxjgj(a;) (the level sets of which are depicted in Figure 121 right) measures the 
maximum distance from x to any of the input points, so a Euclidean ball of radius qs{x) centered 
at x will enclose all the points and is the smallest ball centered at x that encloses all the points. 

If we form a quasiconvex program from the functions qi, the solution to the program consists of 
a pair (A, x) where A = qs{x) and A is as small as possible. That is, the ball with radius A centered 
at x is the smallest enclosing ball of the input points. 

Any smallest enclosing ball problem has a basis of at most d + \ points that determine its value. 
More generally, it will turn out that any quasiconvex program's value is similarly determined by 
a small number of the input functions; this phenomenon will prove central in our ability to apply 
generalized linear programming algorithms to solve quasiconvex programs. 

If we generalize each qi to be the Euclidean distance to a convex set Ki, the resulting quasi- 
convex program finds the smallest sphere that touches or encloses each Ki. In a slightly different 
generalization, if we let qi{x) = d{x,pi) + r,, a sphere centered at x with radius qi{x) or larger will 
contain the sphere centered at pi with radius rj. So, solving the quasiconvex program with this 
family of functions qi will find the smallest enclosing ball of a family of balls [42,70]. 

2.3 Hyperbolic Smallest Enclosing Ball 

Although we have defined quasiconvex programming in terms of Euclidean space M", the defini- 
tion involves only concepts such as convexity that apply equally well to other geometries such 
as hyperbolic space H". Hyperbolic geometry (e.g. see [50]) may be defined in various ways; 
for instance by letting H" consist of the unit vectors of R""*"^ according to the inner product 
{x,y) = Yli<n(^iyi) ~ ^nyn, and defining the distance d{x,y) = cosh""*^ {x,y). Angles, congruence, 
lines, hyperplanes, and other familiar Euclidean concepts can also be defined in a straightforward 
way for hyperbolic space. Hyperbolic geometry satisfies many of the same axioms as Euclidean 
geometry, but not the famous parallel postulate: in the hyperbolic plane H^, given a line £ and a 
point p ^ i, there will be infinitely many lines through p that do not meet i. A hyperbolic convex 
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Fig. 3. Poincare (left) and Klein (rigiit) models of tiie hyperbolic plane. Both models show the 
same hyperbolic arrangement of lines; analogous models exist for any higher dimensional hyperbolic 
space. Figure taken from [8]. 

set K is defined as in Euclidean space to be one in which, for any two points {p, q} C K, all points 
on the line segment connecting p to q also belong to K. Similarly, a quasiconvex function H" i— > M 
is one for which all lower level sets are convex, or equivalently one that is unimodal on any line in 
H". As in the Euclidean case we may define a hyperbolic quasiconvex program to be the problem of 
searching for the point minimizing the pointwise maximum of a collection of hyperbolic quasiconvex 
functions. 

There are several standard ways of representing the points and other geometric objects of 
Hyperbolic space within a Euclidean space, of which the two best known are the Poincare and 
Klein models (Figure 13). In the Poincare model, the points of H" are represented as Euclidean 
points interior to an n-dimensional unit ball or halfspace, and lines of are represented as arcs 
of circles that meet the boundary of this unit ball or halfspace perpendicularly. In this model, the 
hyperbolic angle between two objects in is equal to the Euclidean angle between the models 
of those objects, and hyperbolic circles and spheres are modeled by Euclidean circles and spheres; 
however, hyperbolic distances do not equal distances within the Poincare model, and objects that 
are straight or flat hyperbolically may have curved models. In the Klein model, again, points of 
are represented as Euclidean points interior to an n-dimensional unit ball, but the hyperbolic line 
connecting two points is represented as the restriction to the ball of the Euclidean line connecting 
the models of those points. In this model, angles and distances may be distorted but straightness is 
preserved: a straight or flat hyperbolic object will have a straight or flat model. In particular, since 
the definition of convexity involves only straight line segments, a convex hyperbolic object will have 
a convex Klein model and vice versa. The Poincare and Klein models for a hyperbolic space are 
not uniquely defined, as one may choose any hyperbolic point to be modeled by the center of the 
Euclidean unit ball, and that ball may rotate arbitrarily around its center. 

If we let be a function mapping to a Klein model in M", and if each qi{x) is a hyperbolic 
quasiconvex function, then qi{x) = qi{k^^{x)) is a Euclidean quasiconvex function. More, qi has 
bounded lower levels sets since they are all subsets of the unit ball. Let (A, x) be the solution to the 
Euclidean quasiconvex program defined by the set of functions qi. Then, if x is interior to the unit 
ball defining the Klein model, {\,k~^{x)) is the solution to the hyperbolic quasiconvex program 
defined by the original functions qi. On the other hand, x may be on the boundary of the Klein 
model; if so, x may be viewed as an infinite point of the hyperbolic space, and is the limit of sequence 
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of points within the space with monotonically decreasing values. The latter possibility, of an infinite 
solution to the quasiconvex program, can only occur if some of the hyperbolic quasiconvex functions 
have unbounded lower level sets. Therefore, as Bern and Eppstein [8] noted, hyperbolic quasiconvex 
programs may in general be solved as easily as their Euclidean counterparts. 

As an example, consider the problem of finding the hyperbolic ball of minimum radius containing 
all of a collection of hyperbolic points pi. As in the Euclidean case, we can define qi{x) to be the 
(hyperbolic) distance from x to pi] this function has convex hyperbolic balls as its level sets, so it 
is quasiconvex. And, just as in the Euclidean case, the solution to the quasiconvex program defined 
by the functions qi is the pair (A, x) where the hyperbolic ball of radius A centered at x is the 
minimum enclosing ball of the points pi. 

2.4 Optimal Illumination 

Suppose that we have a room (modeled as a possibly nonconvex three-dimensional polyhedron) and 
wish to place a point source of light in order to light up the whole room as brightly as possible: that 
is, we wish to maximize the minimum illumination received on any point of the room's surface. The 
quasiconvex programs we are studying solve min-max rather than max-min problems, but that is 
easily handled by negating the input functions. 

So, let qi{x) be the negation of the intensity of light received at point i of the room's surface, 
as a function of x, the position of the light source. It is not hard to see that, within any face of 
the polyhedron, the light intensity is least at some vertex of the face, since those are the points 
at maximal distance from the light source and with minimal angle to it. Therefore, we need only 
consider a finite number of possibilities for i: one for each pair (/, v) where / is a face of the 
polyhedron and w is a vertex of /. For each such pair, we can compute qi via a simple formula of 
optics, qi{x) = —u ■ {x — v)/d{x,v)'^, where d is as usual the Euclidean distance, and n is a unit 
vector facing inwards at a perpendicular angle to /. In this formula, one factor u ■ (x — v)/d{x,v) 
accounts for the angle of incidence of light from the source onto the part of face v near vertex 
/, while the other factor l/d(x,v)'^ accounts for the inverse-square rule for falloff of light from a 
point source in three-dimensional space. Note that we can neglect occlusions from other faces in 
this formula, because, if some face is occluded, then at least one other face will be facing away from 
the light source and entirely unilluminated; this unilluminated face will dominate the occluded one 
in our min-max optimization. 

In [3] , as part of a proof of quasiconvexity of a more complex function used for smoothing three- 
dimensional meshes by solid angles, we showed that the function qi defined above is quasiconvex; 
more precisely, we showed that {—qi{x))~^^'^ is a convex function of x by using Mathematica to 
calculate the principal determinants of its Hessian, and by showing from the structure of the 
resulting formulae that these determinants are always nonnegative. Therefore, we can express the 
problem of finding an optimal illumination point as a quasiconvex program. 

2.5 Longest Intersecting Prefix 

This example is due to Chan [14]. Suppose we are given an ordered sequence of convex sets Ki, 
< i < n, that are all subsets of the same compact convex set X C M.'^. We would like to find the 
maximum value i such that Di^iKi ^ 0. That is, we would like to find the longest prefix of the 
input sequence, such that the convex sets in this prefix have a nonempty intersection (Figure E}- 

To represent this as a quasiconvex program, define a nested convex family Ki : Z ^ /^(M'^) for 
each set Ki in the sequence, as follows: 




Ki, a X<-i 

X, otherwise. 
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Fig. 4. Instance of a longest intersecting prefix problem. The longest intersecting prefix is 




Fig. 5. Conversion of convex program into quasiconvex program, by treating each halfspace con- 
straint as a quasiconvex step function. 

The optimal value (A, x) for the quasiconvex program formed by this set of nested convex families 
has X G Ki(A) = Ki for all i < —A, so the prefix of sets with index up to (but not including) —A 
has a nonempty intersection containing x. Since the quasiconvex program solution minimizes A, —A 
is the maximum value with this property. That is, the first —A values of the sequence Ki form its 
longest intersecting prefix. 

More generally, the same technique applies equally well when each of the convex sets Ki has 
an associated value ki, and we must find the maximum value i such that r\ki<iKi ^ 0. The longest 
intersecting prefix problem can be seen as a special case of this problem in which the values ki form 
a permutation of the integers from to n — 1. We will see an instance of this generalized longest 
intersecting prefix problem, in which the values ki are integers with some repeated values, when we 
describe Chan's solution to the Tukey median problem. 

2.6 Linear, Convex, Quasiconvex 

There are many ways of modeling linear programs, but one of the simplest is the following: a linear 
program is the search for a vector x that satisfies all of a set of closed linear inequalities cii - x >bi 
and that, among all such feasible vectors, minimizes a linear objective function f{x) = c - x. The 
vectors x, cii, and c all have the same dimension, which we call the dimension of the linear program. 
We typically use the symbol n to denote the number of inequalities in the linear program. It is often 
useful to generalize such programs somewhat, by keeping the linear constraints but allowing the 
objective function f{x) to be convex instead of linear; such a generalization is known as a convex 
program, and many linear programming algorithms can be adapted to handle the convex case as 
well. 

For instance, consider the following geometric problem, which arises in collision detection algo- 
rithms for maintaining simulations of virtual environments: we are given as input two A;-dimensional 
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convex bodies P and Q, specified as intersections of halfspaces P = OPi and Q = CiQi; we wish to 
find the closest pair of points p, q with p ^ P and g E Q. If we view p, q as forming a 2A:-dimensional 
vector X, then each constraint p G Pj or g S is linear in x, but the objective function d{p, q) is 
nonlinear: evaluating the distance using the Pythagorean formula results in a formula that is the 
square root of a sum of squares of differences of coordinates. We can square the formula to elimi- 
nate the square root, but what remains is a convex quadratic function. Thus, the closest distance 
problem can be expressed as a convex program; similar formulations are also possible when P and 
Q are expressed as convex hulls of their vertex sets [66]. 

These formulations seem somewhat different from our quasiconvex programming framework: in 
the linear and convex programming formulations above, we have a large set of constraints and a 
single objective function, while in quasiconvex programming we have many input functions that 
take a role more analogous to objectives than constraints. Nevertheless, as we now show, any linear 
or convex program can be modeled as a quasiconvex program. Intuitively, the idea is simply to 
treat each halfspace constraint as a quasiconvex step function, and include them together with 
the convex objective functions in the set of quasiconvex functions defining a quasiconvex program 
(Figure ^ . 

Theorem 3. Suppose a convex program is specified by n linear inequalities di-x > bi and a convex 
objective function f{x), and suppose that the solution of this convex program is known to lie within 
a compact convex region K. Then we can find a set ofn + 1 nested convex families Ki{\) such that 
the solution (A, x) of the quasiconvex program formed by these nested convex families is an optimal 
solution to the convex program, with A = f{x). 

Proof. For each inequality di- x >bi form a nested convex family Kj(A) = K {x \ ai ■ x > bi}; that 
is, Ki ignores its argument A and produces a constant compact convex set of the points satisfying 
the ith inequality. Also form a nested convex family Kn = i<'f,K representing the objective function. 

If (A, x) is the optimal solution to the quasiconvex program defined by the nested convex families 
Kj, then di ■ x > bi (else x would not be contained in Ki(A)) and A = f{x) (else either x would 
be outside Kn{X) or the pair {f{x),x) would be a better solution). There could be no y satisfying 
all constraints di ■ y > bi with f{y) < A, else {f{y),y) would be a better solution than (A, x) for 
the quasiconvex program. Therefore, x provides the optimal solution to the convex program as the 
result states. □ 

The region K is needed for this result as a technicality, because our quasiconvex programming 
formulation requires the nested convex families to be compact. In practice, though, it is not generally 
difficult to find K; for instance, in the problem of finding closest distances between convex bodies, we 
could let K he a bounding box defined by extreme points of the convex bodies in each axis-aligned 
direction. 

3 Algorithms 

We now discuss techniques for solving quasiconvex programs, both numerically and combinatorially. 
3.1 Generalized Linear Programming 

Although linear programs can be solved in polynomial time, regardless of dimension [52,53], known 
results in this direction involve time bounds that depend not just on the number and dimension 
of the constraints, but also on the magnitude of the coordinates used to specify the constraints. 
In typical computational geometry applications the dimension is bounded but these magnitudes 
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may not be, so there has been a long hne of work on linear programming algorithms that take 
a linear amount of time in terms of the number of constraints, independent of the magnitude of 
coordinates, but possibly with an exponential dependence on the dimension of the problem [1, IS- 
IS, 28, 66, 68, 69, 83]. In most cases, these algorithms can be interpreted as dual simplex m,ethods: 
as they progress, they maintain a basis of d constraints, and the point x optimizing the objective 
function subject to the constraints in the basis. At each step, the basis is replaced by another one 
with a worse value of x; when no more basis replacement steps are possible, the correct solution 
has been found. 

Very quickly, workers in this area realized that similar techniques could also be applied to certain 
nonlinear programs such as the minimum enclosing ball problem [1,2,15,18,26,27,35,40,41,66,68, 
76,92]. One of the most popular and general formulations of this form of generalized linear program 
is the class of LP-type problems defined by Matousek et al. [66]; we follow the description of this 
formulation from Amenta et al. [3]. 

An LP-type problem consists of a finite set S of constraints and an objective function f mapping 
subsets of S to some totally ordered space and satisfying the following two properties: 

1. For any AcB, f{A) < f{B). 

2. For any A, p, and q, if f{A) = f{A U {p}) = f{A U {q}), then f{A) = f{A U {p, q}). 

The problem is to compute f{S) using only evaluations of / on small subsets of S. 

For instance, in linear programming, S* is a set of halfspaces and f{S) is the point in the 
intersection of the halfspaces at which some linear function takes its minimum value. In the smallest 
enclosing ball problem, S consists of the points themselves, and f(A) is the smallest enclosing ball 
of A, where the total ordering on balls is given by their radii. It is not hard to see that this system 
satisfies the properties above: removing points can only make the radius shrink or stay the same, 
and if a ball contains the additional points p and q separately it contains them both together. 

A basis of an LP-type problem is a set B such that for any AC B, f{A) < f{B). Thus, due to 
the first property of an LP-type problem, the value of the overall problem is the same as the value 
of the optimal basis, the basis B that maximizes f{B). The dimension of an LP-type problem is the 
maximum cardinality of any basis; although we have not included it above, a requirement that this 
dimension be bounded is often included in the definition of an LP-type problem. The dimension 
of an LP-type problem may differ from the dimension of some space that may be associated in 
some way with the problem; for instance, for smallest enclosing balls in R'^, the dimension of the 
LP-type problem turns out to be d + 1 instead of d. 

As Matousek et al. [66] describe, efficient and simple randomized algorithms for bounded- 
dimension LP-type problems are known, with running time 0{dnT + t{d)E\ogn) where n is the 
number of constraints, T measures the time to test whether f{B) = /(i?U{x}) for some basis B and 
element x G S*, t{d) is exponential or subexponential, and E is the time to perform a basis-change 
operation in which we must find the basis of a constant-sized subproblem and use it to replace 
the current basis. It is also possible with certain additional assumptions to solve these problems 
detcrministically in time linear in n [15]. 

As Amenta et al. [3] showed, quasiconvex programs can be expressed as LP-type problems, in 
such a way that the dimension of the LP-type problem is not much more than the dimension of 
the domain of the quasiconvex functions; therefore, quasiconvex programs can be solved in a linear 
number of function evaluations and a sublinear number of basis-change operations. 

In order to specify the LP-type dimension of these problems, we need one additional definition: 
suppose we have a nested convex family If Ki(\) does not depend on A, we say that Ki is constant; 
such constant families arose, for instance, in our treatment of convex programs. Otherwise, suppose 
Ki is associated with a quasiconvex function gj. If there is no open set S such that qi is constant 
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over S, and if K(t') is contained in the interior of K.{t) for any t' < t, we say that k is continuously 
shrinking. We note that this property is different from the related and more well-known property 
of strict quasiconvexity (a quasiconvex function is strictly quasiconvex if, whenever it is constant 
on a line segment, it remains constant along the whole line containing the segment): Li distance 
from the origin (in R'^, d > 1) is continuously shrinking but not strictly quasiconvex. On the other 
hand, the function 

f{x, y) = min{r | + (y — r)^ < r^} 

(on the closed upper halfplane y > 0) is strictly quasiconvex but not continuously shrinking, since 
the origin is on the boundary of all its level sets. 

We repeat the analysis of Amenta et al. [3], showing that quasiconvex programs are LP-type 
problems, below. 

Theorem 4. Any quasiconvex program forms an LP-type problem of dimension at most 2d + l. If 
each Ki in the quasiconvex program is either constant or continuously shrinking, the dimension is 
at most d+ 1. 

Proof. We form an LP-type problem in which the set S consists of the nested convex families defin- 
ing the quasiconvex program, and the objective function A(T) gives the value of the quasiconvex 
program defined by the nested convex families in T. Then, property 1 of LP-type problems is ob- 
vious: adding another nested convex family to the input can only further constrain the solution 
values and increase the min-max solution. To prove property 2, recall that A{T) is defined as the 
minimum point of the intersection {(A, x) \ x € Kj(A)} (the intersection is nonempty by the remark 
in Section [1.31 about replacing infima by minima). If this point belongs to the intersection for sets 
A, Au {kj}, and A U {/tfe}, then clearly it belongs to the intersection for A U {ni, Kj}. It remains 
only to show the stated bounds on the dimension. 

First we prove the dimension bound for the general case, where we do not assume continuous 
shrinking of the families in S. Let (A,x) = A{S). For any A' < A, 

n ^d>^') = 0' 

SO by Helly's theorem some {d + l)-tuple of sets Kj(A') has empty intersection. If there is some 
A" < A for which this {d+ l)-tuple's intersection becomes nonempty, replace A' by A", find another 
{d + l)-tuple with empty intersection for the new A', and repeat until this replacement process 
terminates. There are only finitely many possible {d + l)-tuples of nested convex families, and 
each replacement increases A', so the replacement process must terminate and we eventually find a 
{d + l)-tuple B~ of nested convex families that has empty intersection for all A' < A. 

With this choice of B~ , A{B^) = (A,y) for some y, so the presence of forces the LP-type 
problem's solution to have the correct value of A. We must now add further nested convex families 
to our basis to force the solution to also have the correct value of x. Recall that 

X e L = Ki{X), 

and X is the minimal point in L. By Helly's theorem again, the location of this minimal point is 
determined by some d-tuple B^ of the sets Kj(A). Then A{B^ U B~^) = A(S), so some basis of S is 
a subset of B~ U B~^ and has cardinality at most 2d + 1. 

Finally, we must prove the improved dimension bound for well-behaved nested convex families, 
so suppose each Hi £ S is constant or continuously shrinking. Our strategy will be to again find a 
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tuple B~ that determines A, and a tuple that determines x, but we will use continuity to make 
the sizes of these two tuples add to at most d+1. 

The set L defined above has empty interior: otherwise, we could find an open region X within 
L, and a nested family Hi £ S such that Ki(A') H X = for any A' < A, violating the assumption 
that Ki is constant or continuously shrinking. If the interior of some Ki{X) contains a point of the 
affine hull of L, we say that Ki is "slack" ; otherwise we say that Ki is "tight" . The boundary of a 
slack Ki{\) intersects L in a subset of measure zero (relative to the affine hull of L), so we can find 
a point y in the relative interior of L and not on the boundary of any slack Ki. Form the projection 
TT : M"* 1-^ j^d-dimL QYito the orthogonal complement of L. 

For any ray r in M'^-dimi starting at the point 7r(L), we can lift that ray to a ray f in M*^ 
starting at y, and find a hyperplane containing L and separating the interior of some Ki{X) from 
f \ {y}. This separated Kj must be tight (because it has y on its boundary as the origin of the ray) 
so the separating hyperplane must contain the affine hull of L (otherwise some point in L within a 
small neighborhood of x would be interior to Ki). Therefore the hyperplane is projected by tt to a 
lower dimensional hyperplane separating 7r(/«i(A)) from 7r(L). Since one can find such a separation 
for any ray, HtightK; 7r(Kj(A)) can not contain any points of any such ray and must consist of the 
single point 7r(L). At least one tight Kj must be continuously shrinking (rather than constant), 
since otherwise f]^.^s i^ii^') would be nonempty for some A' < A. The intersection of the interior of 
Tr{Kj{\)) with the remaining projected tight constraints Tr{Ki{X)) is empty, so by Helly's theorem, 
wc can find a (d — dimL + l)-tuple B~ of these convex sets having empty intersection, and the 
presence of forces the LP-type problem's solution to have the correct value of A. Similarly, we 
can reduce the size of the set Z?+ determining x from d to dimL, so the total size of a basis is at 
most (d — dim L + 1) + dim L = d+1. □ 

This result provides theoretically efficient combinatorial algorithms for quasiconvex programs, 
and allows us to claim 0{n) time randomized algorithms for most quasiconvex programming prob- 
lems in the standard computational model for computational geometry, in which primitives of 
constant description complexity may be assumed to be solved in constant time. For certain well- 
behaved sets of quasiconvex functions (essentially, the family of sets Sx,\ = {k £ S \ x € «^(A)} 
should have bounded Vapnik-Chervonenkis dimension) the technique of Chazelle and Matousek [15] 
applies and these problems can be solved deterministically in 0{n) time. 

However, we should note that there are some difficulties with this approach in practice. In partic- 
ular, although the basis-change operations have constant description complexity, it may not always 
be clear how to implement them efficiently. Therefore, in the next section we discuss alternative 
numerical techniques for solving quasiconvex programs directly, based only on simpler operations 
(function and gradient evaluation). It may be of interest to combine the two approaches, by us- 
ing numerical techniques to solve the basis change operations needed for the LP-type approach; 
however, we do not have any theory describing how the LP-type algorithms might be affected by 
approximate numerical results in the basis-change steps. 

3.2 Implicit Quasiconvex Programming 

In some circumstances we may have a set of n inputs that leads to a quasiconvex program with 
many more than n quasiconvex functions; for instance, there may be one such function per pair of 
inputs. If we directly apply an LP-type algorithm, we will end up with a running time much larger 
than the 0{n) input size. Chan [14] showed that, in such circumstances, the time for solving the 
quasiconvex program can often be sped up to match the time for a decision algorithm that merely 
tests whether a given pair (A, x) provides a feasible solution to the program. 
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As a simple example, consider a variation of the smallest enclosing ball problem. Suppose that 
we wish to place a center that minimizes the maximum sum of distances to any fe-tuple of sites, 
rather than (as in the smallest enclosing ball problem) minimizing the maximum distance to a 
single site. This can be expressed again as a quasiconvex program: the sum of distances to any 
fe-tuple of sites is quasiconvex, as it is a sum of convex functions. There are 0{n^) such functions, 
so the problem can be solved in 0{n^) time by the methods discussed already. However, the quality 
of any fixed center can easily be evaluated much more quickly, in 0{n) time, and Chan's technique 
provides an automatic method for turning this fast evaluation algorithm into a fast optimization 
algorithm for choosing the best center location. 

Chan's result applies more generally to LP-type problems, but we state it here as it applies to 
implicit quasiconvex programming. 

Theorem 5. Let Q be a space of quasiconvex functions, V be a space of input values, and / : 2^ i— > 
map sets of input values to sets of functions in Q. Further, suppose that V, f, and S satisfy 
the following properties: 

— There exists a constant-time subroutine for solving quasiconvex programs of the form f{B) for 
any B dV with \B\ = 0(1). 

— There exists a decision algorithm that takes as input a set P dV and a pair (A, x), and returns 
yes if and only if x E k{X) for all k G f{P)- The running time of the decision algorithm 
is hounded by D{\P\), where there exists a constant e > such that D{n)/n^ is monotone 
increasing. 

— There are constants a and r such that, for any input set P C V, we can find in time at most 
D{\P\) a collection of sets Pi, < i < r, each of size at most a\P\, for which f{P) = Uif{Pi). 

Then for any P C V we can solve the quasiconvex program f{P), where \P\ = n, in randomized 
expected time 0{D{n)). 

The proof involves solving a slightly more general problem in which we are given, not just a 

single input P, but a set of inputs Pi, . . ., P^, where d is the dimension of the LP-type problems 
coming from Q, and must solve the quasiconvex program U/(Pj). Given any such problem, we 
partition each input Pi into r* subproblems Pij of size at most Q*n for an appropriately chosen 
i, by repeatedly subdividing large subproblems into smaller ones. We then view the subproblems 
Pij as being constraints for an LP-type problem in which the objective function is the solution 
to the quasiconvex program Up. .^sfiPi,])- This new LP-type problem turns out to have the same 
dimension as the quasiconvex programs with which we started, and the result follows by applying 
a standard LP-type algorithm to this problem and solving the divide-and-conquer recurrence that 
results. 

The first and last conditions of the theorem are easily met when f{P) produces one or a constant 
number of quasiconvex functions per fe-tuple of inputs for some constant k (as in our example of 
optimizing the sum of k distances): then, constant sized input sets lead to constant sized quasiconvex 
programs, and if the input is partitioned into k + 1 equal-sized subsets, the complements of these 
subsets provide the sets Pi needed for the last condition. For such problems, the main difficulty in 
applying this theorem is finding an appropriate decision algorithm. For our example of minimizing 
the maximum sum of k distances, the decision algorithm is also straightforward (select and add the 
k largest distances from the given center to the sites) and so we can apply Chan's result to solve 
this problem in 0{n) time. 

Chan's implicit quasiconvex programming algorithm is important in the robust statistics ap- 
plication described later. This algorithm has also been applied to problems of inverse parametric 
minimum spanning tree computation [14,31] and facility location [34]. 
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Fig. 6. Example showing the difficulty of applying standard gradient descent methods to quasi- 
convex programming. The function to be minimized is the maximum distance to any point; only 
points within the narrow shaded intersection of circles have function values smaller than the value 
at point w. Figure taken from [33]. 



3.3 Smooth Quasiconvex Programming 

If all functions qi{x) are quasiconvex, the function q{x) = maxj qi{x) is itself quasiconvex, so we can 
apply hill-climbing procedures to find its infimum. Such hill climbing procedures may be desirable 
in preference to the combinatorial algorithms for LP-type problems, as they avoid the difficulty of 
describing and implementing an appropriate exact basis change procedure. In addition, a hill climb- 
ing information that uses only numerical evaluation of function values (or possibly also function 
gradient evaluations) can be implemented in a generic way that does not depend on the specific 
form of the quasiconvex functions given to it as input. 

However, many of the known non-linear optimization techniques require the function being op- 
timized to satisfy some smoothness conditions. In many of our applications the individual functions 
Qi are smooth, but their maximum q may not be smooth, so it is difficult to apply standard gradient 
descent techniques. The difficulty may be seen, for instance, in the smallest enclosing ball problem 
in the plane (Figure IHl). A basis for this problem may consist of either two or three points. If a 
point set has only two points in its basis, and our hill climbing procedure for circumradius has 
reached a point w equidistant from these two points and near but not on their midpoint, then 
improvements to the function value q{w) may be found only by moving w in a narrow range of 
directions towards the midpoint. Standard gradient descent algorithms may have a difficult time 
finding such an improvement direction. 

To avoid these difficulties, we introduced in [33] the following algorithm, which we call smooth 
quasiconvex programming, and which can be viewed as a generalization of Zoutendijk's method of 
feasible directions [94] for convex programming. If a quasiconvex function qi is differentiable, and 
t(; is a point where qi is not minimal, then one can find a point with a smaller value by moving a 
sufficiently small distance from x along any direction having negative dot product with the gradient 
of qi at w. Thus, we can improve q{w) by moving in a direction that is negative with respect to all 
the gradients of the functions that determine the value of q{w). 

We formalize this notion and generalize it to nondifferentiable functions as follows. Assume for 
the purposes of this algorithm that, for each of the input quasiconvex functions qi, and each x that 
is not the minimum point of qi, we also can compute a vector- valued function g*(a;), satisfying the 
following properties: 

1. If qi{y) < qi{x), then {y - x) ■ q*{x) > 0, and 
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2. If q*{x) • y > 0, then for all sufficiently small e > 0, qi{x + ey) < qi{x). 

Less formally, any vector y is an improving direction for qi{x) if and only if it has positive inner 
product with ql{x)- 

If the level set qf^ is a smooth convex set (one that has at each of its boundary points a unique 
tangent plane), then the vector q*{x) should be an inward-pointing normal vector to the tangent 
plane to q^ at x. For example, in the smallest enclosing ball problem, the level sets arc spheres, 
having tangent planes perpendicular to the radii, and q* should point inwards along the radii of 
these spheres. If qi is differentiable then q^ can be computed as the negation of the gradient of qi, 
but the functions also exist for discontinuous functions with smooth level sets. 

Our smooth quasiconvex programming algorithm begins by selecting an initial value for x, and 
a desired output tolerance. Once these values are selected, we repeat the following steps: 

1. Compute the set of vectors q*{x), for each i such that qi{x) is within the desired tolerance of 
maxigj(x). 

2. Find an improving direction y; that is, a vector such that y ■ q^(x) > for each vector q^{x) in 
the computed set. If no such vector exists, q{x) is within the tolerance of its optimal value and 

the algorithm terminates. 

3. Search for a value e for which q{x + ey) < q{w), and replace x hy x + ey. 

The search for a vector y in step 2 can be expressed as a linear program. However, when the 
dimension of the quasiconvex functions' domain is at most two (as in the planar smallest enclosing 
ball problem) it can be solved more simply by sorting the vectors q*{x) radially around the origin 
and choosing y to be the average of two extreme vectors. 

In step 3, it is important to choose e carefully. It would be natural, for instance, to choose e 
as large as possible while satisfying the inequality in that step; such a value could be found by a 
simple doubling search. However, such a choice could lead to situations where the position of x 
oscillates back and forth across the true optimal location. Instead, it may be appropriate to reduce 
the resulting e by a factor of two before replacing x. 

We do not have any theory regarding the convergence rate of the smooth quasiconvex program- 
ming algorithm, but we implemented it and applied it successfully in the automated algorithm 
analysis application discussed below [33]. Our implementation appeared to exhibit linear conver- 
gence: each iteration increased the number of bits of precision of the solution by a constant. Among 
numerical algorithms techniques, the sort of gradient descent we perform here is considered naive 
and inefficient compared to other techniques such as conjugate gradients or Newton iteration, and 
it would be of interest to see how well these more sophisticated methods could be applied to 
quasiconvex programming. 

4 Applications 

We have already described some simple instances of geometric optimization problems that can be 
formulated as quasiconvex programs. Here wc describe some more complex applications of geometric 
optimization, in which quasiconvex programming plays a key role. 

4.1 Mesh Smoothing 

An important step in many scientific computation problems, in which differential equations de- 
scribing airflow, heat transport, stress, global illumination, or similar quantities are simulated, is 
mesh generation [7,10]. In this step, a complex two- or three-dimensional domain is partitioned into 
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Fig. 7. Mesh of an arched domain. Too much Laplacian smoothing can lead to invahd placements 
of the internal vertices beyond the boundaries of the arch. 




Fig. 8. Optimization-based smoothing of a triangular mesh in M^. At each step we remove a vertex 
from the mesh, leaving a star-shaped polygon, then add a new vertex within the kernel (shaded) 
of the star-shaped region and retriangulate. Figure taken from [3]. 

simpler regions, called elements, such as triangles or quadrilaterals in the plane or tetrahedra or 
cuboids in three dimensions. Once these elements are formed, one can then set up simple equations 
relating the values of the quantity of interest in each of the elements, and solve the equations to 
produce the results of the simulation. In this section we are particularly concerned with unstruc- 
tured mesh generation, in which the pattern of connections from element to element does not form 
a regular grid; we will consider a problem in structured mesh generation in a later section. 

In meshing problems, it is important to find a mesh that has small elements in regions of fine 
detail, but larger elements elsewhere, so that the total number of elements is minimized; this allows 
the system of equations derived from the mesh to be solved quickly. It is also important for the 
accuracy of the simulation that the mesh elements be well shaped; typically this means that no 
element should have very sharp angles or angles very close to 180°. To achieve a high quality mesh, 
it is important not only to find a good initial placement of mesh vertices (the main focus of most 
meshing papers) but then to modify the mesh by changing its topology and moving vertices until 
no further quality increase can be achieved. We here concentrate on the problem of moving mesh 
vertices while retaining a fixed mesh topology, known as mesh smoothing [3,4,13,24,36-39,91]. 

Two approaches to mesh smoothing have commonly been used, although they may sometimes 
be combined [13,36]: In Laplacian smoothing, all vertices are moved towards the centroid of their 
neighbors. Although this is easy and works well for many instances, it has some problems; for 
instance in a regular mesh on an arched domain (Figure E|) , repeated Laplacian smoothing can 
cause the vertices at the top of the arch to sag downwards, eventually moving them to invalid 
positions beyond the boundaries of the domain. 

Instead, optimization-based smoothing takes a more principled approach, in which we decide on 
a measure of element quality that best fits our application, and then seek the vertex placement 
that optimizes that quality measure. However, since simultaneous global optimization of all vertex 
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min max area 
min max altitude 



max min area 
max min altitude 
min max aspect ratio 



min max angle 



max min angle max min angle min max aspect ratio 

max min altitude 



min max perimeter max min edge length min max enclosing disk 
min max diameter 

Fig. 9. Level set shapes for various mesh element quality measures. Figure modified from one in [3]. 



positions seems a very difficult problem, we instead cycle through the vertices optimizing their 
positions a single vertex at a time. At each step (Figure |H1), we select a vertex and remove it from 
the mesh, leaving a star-shaped region consisting of the elements incident to that vertex. Then, 
we place a new vertex within the kernel of the star-shaped region, and form a mesh again by 
connecting the new vertex to the boundary of the region. Each step improves the overall mesh 
quality, so this optimization process eventually converges to a locally optimal placement, but we 
have no guarantees about its quality with respect to the globally optimal placement. 

However, in the individual vertex placement steps we need accept no such compromises with 
respect to global optimization. As we showed in [3], for many natural measures qi{x) of the quality 
of an element incident to vertex x (with smaller numbers indicating better quality), the problem of 
finding a mesh minimizing the maximum value of qi can be expressed as a quasiconvex program. Fig- 
ured illustrates the level set shapes resulting from various of these quasiconvex optimization-based 
mesh smoothing problems. For shape-based quality measures, such as maximizing the minimum 
angle, the optimal vertex placement will naturally land in the interior of the kernel of the region 
formed by the removal of the previous vertex placement. For some other quality measures, such as 
minimizing the maximum perimeter, it may be appropriate to also include constant quasiconvex 
functions, forcing the vertex to stay within the kernel, similar to the functions used in our transfor- 
mation of convex programs to quasiconvex programs. It would also be possible to handle multiple 
quality measures simultaneously by including quasiconvex functions of more than one type in the 
optimization problem. 

In most of the cases illustrated in Figure El it is straightforward to verify that the quality 
measure has level sets of the convex shape illustrated. One possible exception is the problem of 
minimizing the maximum aspect ratio (ratio of the longest side length to shortest altitude) of any 
element. To see that this forms a quasiconvex optimization problem. Amenta et al. [3] consider 
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separately the ratios of the three sides to their corresponding altitudes; the maximum of these 
three will give the overall aspect ratio. The ratio of a side external to the star to its corresponding 
altitude has a feasible region (after taking into account the kernel constraints) forming a halfspace 
parallel to the external side, as shown in Figure IH] (top center). To determine the aspect ratio on 
one of the other two sides of a triangle Ai, normalize the triangle coordinates so that the replaced 
point has coordinates (x, y) and the other two have coordinates (0,0) and (1,0). The side length is 
then Y^x^ + y^, and the altitude is yl\Jx^ + y^, so the overall aspect ratio has the simple formula 
(x^ +y^) jy. The locus of points for which this is a constant h is given by = 6y, or equivalently 

+ (y ~ ^I'^yf' = (^/2)^- Thus the feasible region is a circle tangent to the fixed side of Ai at 
one of its two endpoints (Figure IHl center right). Another nontrivial case is that of minimizing the 
smallest enclosing ball of the element, shown in the bottom right of the figure; in that case the level 
set boundary consists of curves of two types, according to whether, for placements in that part of 
the level set, the enclosing ball touches two or three of the element vertices, but the curves meet 
at a common tangent point to form a smooth convex level set. 

Bank and Smith [4] define yet another measure of the quality of a triangle, computed by dividing 
the triangle's area by the sum of the squares of its edge lengths. This gives a dimensionless quantity 
which Bank and Smith normalize to be one for the equilateral triangle (and less than one for any 
other triangle). As Bank and Smith show, the lower level sets for this mesh quality measure form 
circles centered on the perpendicular bisector of the two fixed points of the mesh element, so, as with 
the other measures, finding the placement optimizing Bank and Smith's measure can be expressed 
as a quasiconvex program. 

We have primarily discussed triangular mesh smoothing here, but the same techniques apply 
with little modification to many natural element quality measures for quadrilateral and tetrahedral 
mesh smoothing. Smoothing of cubical meshes is more problematic, though, as moving a single 
vertex may cause the faces of one of the cuboid elements to become significantly warped. Several 
individual quasiconvex quality measures for quadrilateral and tetrahedral meshes, and the shapes 
of their level sets, are discussed in more detail in [3]. The most interesting of these from the 
mathematical viewpoint is the problem of maximizing the minimum solid angle of any tetrahedral 
element, as measured at its vertices, which with some difficulty we were able to show leads to a 
quasiconvex objective function. 

4.2 Graph Dravi^ing 

The Koebe-Thurston-Andreev embedding theorem [11,55,81] states that any planar graph embed- 
ding can be transformed into a collection of disks with disjoint interiors on the surface of a sphere, 
one disk per vertex, such that two disks are tangent if and only if the corresponding two vertices 
are adjacent (Figure [Till left and center). The representation of the graph as such a collection of 
tangent disks is sometimes called a coin graph. For maximal planar graphs, this coin graph repre- 
sentation is unique up to Mobius transformations (the family of transformations of the sphere that 
transform circles to circles), and for non- maximal graphs it can be made unique by adding a new 
vertex within each face of the embedding, adjacent to all vertices of the face, and finding a disk 
representation of the resulting augmented maximal planar graph. 

Given a coin graph representation, the graph itself can be drawn on the sphere e.g. by placing 
a vertex at the center of each circle and connecting two vertices by edges along an arc of a great 
circle; similar drawings are also possible in the plane by using polar projection to map the circles in 
the sphere onto circles in the plane [45]. Coin graphs can also be used to form a three-dimensional 
polyhedral representation of the graph, as follows: embed the sphere in space, and, for each disk, 
form a cone in space that is tangent to the sphere at the disk's boundary; then, form a polyhedron 
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Fig. 10. Planar graph (left), its representation as a set of tangent disks on a sphere (center), and 
the corresponding polyhedral representation (right). Left and center images taken from [8]. 



by taking the convex hull of the cone apexes. The resulting polyhedron's skeleton is isomorphic to 
the original graph, and its edges are tangent to the sphere (Figure ITUl right). 

In order to use these techniques for visualizing graphs, we would like to choose a coin graph 
representation that leads to several desirable properties identified as standard within the graph 
drawing literature [5], including the display of as many as possible of the symmetries of the original 
graph, and the separation of vertices as far apart from each other as possible. Our paper with 
Bern [8] used quasiconvex programming to formalize the search for a drawing based on these 
objectives. 

In order to understand this formalization, we need some more background knowledge about 
Mobius transformations and their relation to hyperbolic geometry. We can identify the unit sphere 
that the Mobius transformations transform as being the boundary of a Poincare or Klein model of 
hyperbolic space H^. The points on the sphere can be viewed as "infinite" points that do not belong 
to but are the limit points of certain sequences of points within H^. With this identification, 
circles on the sphere become the limit points of hyperplanes in H^. Any isometry of takes 
hyperplanes to hyperplanes, and therefore can be extended to a transformation of the sphere that 
takes circles to circles, and the converse turns out to be true as well. We can determine an isometry 
of by specifying which point of is mapped to the center of the Poincare or Klein model, 
and then by specifying a spatial rotation around that center point. The rotation component of this 
isometry does not change the shape of objects on the sphere, so whenever we seek the Mobius 
transformation that optimizes some quality measure of a transformed configuration of disks on the 
sphere, we can view the problem more simply as one of seeking the optimal center point of the 
corresponding isometry in H'^. 

To see how we apply this technique to our graph drawing problem, first consider a version of 
the problem in which we seek a disk representation maximizing the radius of the smallest disk. 
More generally, given any collection of circles on the sphere, we wish to transform the circles in 
order to maximize the minimum radius. Thus, let qi{x) measure the (negation of the) transformed 
radius of the ith circle, as a function of the transformed center point x G H'^. If we let Hi denote 
the hyperplane in that has the ith circle as its set of limit points, then the transformed radius 
is maximized when the circle is transformed into a great circle; that is, when x E Hi. If we choose 
a center point x away from Hi, the transformed radius will be smaller, and due to the uniform 
nature of hyperbolic space the radius can be written as a function only of the distance from x 
to Hi, not depending in any other way on the location of x. That is, the level sets of qi are the 
convex hyperbolic sets within some distance R of the hyperplane H^. Therefore, qi is a quasiconvex 
hyperbolic function. In fact, the quasiconvex program defined by the functions qi can be viewed as 
a hyperbolic version of a generalized minimum enclosing ball problem, in which we seek the center 
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Fig. 11. Two-dimensional analogue of max-min radius transform problem: find the smallest disk 
touching all of a collection of hyperbolic lines. 

X of the smallest ball that touches each of the convex sets Hi. The two-dimensional version of this 
problem, in which we seek the smallest disk touching each of a collection of hyperbolic lines, is 
illustrated in Figure 1111 If we form a Klein or Poincare model with the resulting optimal point 
X at the center of the model, the corresponding Mobius transformation of the model's boundary 
maximizes the minimum radius of our collection of circles. 

Further, due to the uniqueness of quasiconvex program optima, the resulting disk representation 
must display all the symmetries possible for the original planar graph embedding; for, if not all 
symmetries were displayed, one could use an undisplayed symmetry to relabel the vertices of the 
disk representation, achieving a second disk representation with equal quality to the first. For 
instance, in Figure IIUI the disk representation shown has three planes of mirror symmetry while 
the initial drawing has only one mirror symmetry axis. 

Bern and Eppstein [8] then consider an alternative version of the graph drawing problem, in 
which the objective is to maximize the minimum distance between certain pairs of vertices on 
the sphere surface. For instance, one could consider only pairs of vertices that are adjacent in the 
graph, or instead consider all pairs; in the latter case we can reduce the number of pairs that need 
be examined by the algorithm by using the Delaunay triangulation in place of the complete graph. 
The problem of maximizing the minimum spherical distance among a set of pairs of vertices can 
be formulated as a quasiconvex program by viewing each pair of vertices as the two limit points of 
a hyperbolic line in M^, finding the center x of the smallest ball in that touches each of these 
hyperbolic lines, and using this choice of center point to transform the sphere. 

Mobius transformations can also be performed on the augmented plane U {oo} instead of on 
a sphere, and act on lines and circles within that plane; a line can be viewed as a limiting case of 
a circle that passes through the special point oo. Multiplication of each coordinate of each point 
by the same constant k forms a special type of Mobius transformation, which (if A; > 1) increases 
every distance, so it does not make sense to look for an unrestricted Mobius transformation of 
the plane that maximizes the minimum Euclidean distance among a collection of pairs of points. 
However, Bern and Eppstein were able to show, given a collection of points within the unit ball in 
the plane, that seeking the Mobius transformation that takes that disk to itself and maximizes the 
minimum distances between certain pairs of the points can again be expressed as a two-dimensional 
quasiconvex program. The proof of quasiconvexity is more complex and involves simultaneously 
treating the unit ball as a Poincare model of and the entire plane as the boundary of a Poincare 
model of H^. 
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Fig. 12. Conformal meshing: transform domain to a more simply shaped region with a known mesh, 
then invert the transformation to transform the mesh back to the original domain. 

Along with these coin graph based drawing methods, Bern and Eppstein also considered a 
different graph drawing question, more directly involving hyperbolic geometry. The Poincare and 
Klein models of projective geometry have been considered by several authors as a way of achieving 
a "fisheye" view of a large graph, so that a local neighborhood in the graph is visible in detail 
near the center of the view while the whole graph is spread out on a much smaller scale at the 
periphery [56,71,72]. Bern and Eppstein [8] found quasiconvex programming formulations of several 
versions of the problem of selecting an initial viewpoint for these hyperbolic drawings, in order for 
the whole graph to be visible in as large a scale as possible. For instance, a natural version of 
this problem would be to choose a viewpoint minimizing the maximum hyperbolic distance to any 
vertex, which is just the hyperbolic smallest enclosing ball problem again. One question in this area 
that they left open is whether one can use quasiconvex programming to find a Klein model of a 
given graph that maximizes the minimum Euclidean distance between adjacent vertices. 

4.3 Conformal Mesh Generation 

The ideas of mesh generation and optimal Mobius transformation coincide in the problem of con- 
formal mesh generation [8]. In this problem, we wish to generate a mesh for a simply-connected 
domain in by using a conformal transformation (that is, a transformation that preserves angles 
of incidence between transformed curves) to map the shape into some easy-to-mesh domain such as 
a square, then invert the transformation to map the meshed square back into the original domain 
(Figure [T2|) . There has been much work on algorithms for finding conformal maps [25,47,84,85,89] 
and conformal meshes have significant advantages: the orthogonality of the angles at mesh ver- 
tices means that one can avoid certain additional terms in the definition of the partial differential 
equation to be solved [10,88]. 

If we replace the square in Figure El by a disk, the Riemann mapping theorem tells us that 
a conformal transformation always exists and is, moreover, unique up to Mobius transformations 
that transform the disk to itself; any such transformation preserves conformality. Thus, we have 
several degrees of freedom for controlling the size of the mesh elements produced by the conformal 
method: we can use a larger or smaller grid on the disk or square, but we can also use a Mobius 
transformation in order to enlarge certain portions of the domain and shrink others before meshing 
it. We would like to use these degrees of freedom to construct a mesh that has small elements in 
regions of the domain where fine detail is desired, and large elements elsewhere, in order to limit 
the total number of elements of the resulting mesh. 

Bern and Eppstein [8] formalized the problem by assuming an input domain in which certain 
interior points pi are marked with a desired element size Sj. If we find a conformal map / from 
the domain to a disk, the gradient of / maps the marked element sizes to desired sizes s[ in the 
transformed disk: s[ = \\f'{pi)\\- We can then choose a structured mesh with element size mins^ 
in the disk, and transform it back to a mesh of the original domain. The goal is to choose our 
conformal map in a way that maximizes mins^, so that we can use a structured mesh with as 
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few elements as possible. Another way of interpreting this is that s[ can be seen as the radius of 
a small disk at f{pi). What we seek is the transformation that maximizes the minimum of these 
radii. This is not quite the same as the max-min radius graph drawing problem of the previous 
section, because the circles to be optimized belong to instead of to a sphere, but as in the 
previous section we can view the unit disk as being a Poincare model of (using the fact that 
circles in are mapped by the Poincare model into circles in the unit disk), and seek a hyperbolic 
isometry that maps into itself and optimizes the circle radii. The transformed radius of a circle 
is a function only of the distance from that circle to the center point of the transformed model, 
so the level sets of the functions representing the transformed radii are themselves circles and the 
functions are quasiconvex. 

The quasiconvex conformal meshing technique of Bern and Eppstein does not account for two 
remaining degrees of freedom: first, it is possible to rotate the unit disk around its center point and, 
while that will not change the element size as measured by Bern and Eppstcin's formalization, it 
will change the element orientations. This is more important if we also consider the second degree 
of freedom, which is that instead of using a uniform grid on a square, we could use a rectangle 
with arbitrary aspect ratio. Bern and Eppstein leave as an open question whether we can efficiently 
compute the optimal choice of conformal map to a high-aspect-ratio rectangle to maximize the 
minimum desired element size. 

4.4 Brain Flat Mapping 

Hurdal et al. [49] describe methods for visualizing the complicated structure of the brain by stretch- 
ing its surface onto a flat plane. They perform this stretching via conformal maps: surfaces of major 
brain components such as the cerebellum are simply connected, so there exists a conformal map from 
these surfaces onto a Euclidean unit disk, sphere, or hyperbolic plane. Hurdal et al. approximate 
this conformal map by using a fine triangular mesh to represent the brain surface, and forming the 
Kocbe disk representation of this mesh. Each triangle from the brain surface can then be mapped to 
the triangle connecting the corresponding three disk centers. As in the conformal meshing example, 
there is freedom to modify the conformal map by means of a Mobius transformation, so Bern and 
Eppstein [8] suggested that the optimal Mobius transformation technique described in the previous 
two sections could also be useful in this application. 

Although conformal transformation preserves angles, it distorts other important geometric in- 
formation such as area. Bern and Eppstein proposed to ameliorate this distortion by using an 
optimal Mobius transformation to find the conformal transformation minimizing the maximum ra- 
tio a /a' where a is the area of a triangle in the initial three-dimensional map, and a' is the area of 
its image in the flat map. 

Unfortunately it has not yet been possible to prove that this optimization problem leads to 
quasiconvex optimization problems. Bern and Eppstein formalized the difficulty in the following 
open question: Let T be a triangle in the unit disk or on the surface of a sphere, and let C be 
the set of center points for Poincare models (of in the disk case or H'^ in the sphere case) such 
that the Mobius transformations corresponding to center points in C transform T into a triangle 
of area at least A. Is C necessarily convex? Note that, at least in the spherical case, the area of the 
transformed triangle is the same as the hyperbolic solid angle of T as viewed from the center point, 
so this question seems strongly reminiscent of the difficult problem of proving quasiconvexity for 
tetrahedral mesh smoothing to maximize the minimum Euclidean solid angle, discussed in the initial 
subsection of this section. A positive answer would allow the quasiconvex programming technique 
to be applied to this brain flat mapping application. 
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Fig. 13. An additive color gamut, with vertices labeled by colors: K = black, R = red, G = green, 
B = blue, C = cyan, M = magenta, Y = yellow, W = white. 

4.5 Optimized Color Gamuts 

Tiled projector systems [48, 60, 77] are a recent development in computer display technology, in 
which the outputs of multiple projectors are combined into large seamless displays for collaborative 
workspaces. There are many difficult research issues involved in achieving this seamlessness: how 
to move the data quickly enough to all the screens, how to maintain physical alignment of the 
projectors, how to handle the radial reduction in brightness (vignetting) common to many projector 
systems, and so on. Here we concentrate on one small piece of this puzzle: matching colors among 
the outputs of multiple projectors. Any imaging device has a gamut, the set of colors that it can 
produce. However, two projectors, even of the same model, will have somewhat different gamuts 
due to factors such as color filter batches and light bulb ages. We seek a common gamut of colors 
that can be produced by all the projectors, and a coordinate system for that gamut so that we can 
display color images in a seamless fashion across multiple projectors [9,64,86]. 

Most projectors, and most computer graphics software, use an additive color system in which 
colors are produced by adding signals of three primary colors, typically red, green, and blue. If 
we view the gamuts as sets of points in a linear three-dimensional device-independent color space, 
additive color systems produce gamuts that are the Minkowski sums of three line segments, one per 
color signal, and therefore have the geometric form of parallelepipeds (Figure [TH|) . The color spaces 
representing human vision are three-dimensional, so these parallelepipeds have twelve degrees of 
freedom: three for the black point of the projector (representing the color of light it projects when it 
is given a zero input signal) and three each for the three primary colors (that is, the color that the 
projector produces when given an input signal with full strength in one primary color channel and 
zero in the other two color channels). The black point and the three primary colors form four of 
the eight parallelepiped vertices; the other four are the secondary colors cyan, yellow, and magenta, 
and the white point produced when all three input color channels are saturated. 

The computational task of finding a common color gamut, then, can be represented as a twelve- 
dimensional geometric optimization problem in which we seek the best parallelepiped to use as our 
gamut, according to some measure of gamut quality, while constraining our output parallelepiped 
to lie within the intersection of a collection of input parallelepipeds, one per projector of our system. 

To represent this problem as a quasiconvex program, Bern and Eppstein [9] suppose that we are 
given eight quasiconvex functions dx, dji, dc, ds, dc, djyf, dy, and dw, where each : 
measures the distance of a color from the ideal location of corner X of the color cube (here each 
capital letter is the initial of one of the colors at the color cube corners, except for K which 
by convention stands for black). This formulation allows different distance functions to be used 
for each color; for instance, we might want to weight dx and dw more strongly than the other 
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six color distances. We also form eight functions fx '■ 1^^^ '-^ mapping our twelve-dimensional 
parametrization of color gamuts into the color values of each of the gamut corners. If we parametrize 
a gamut by the black point and three primary colors, then fx, fn, /g) and Jb are simply coordinate 
projections, while the other four functions are simple linear combinations of the coordinates. For 
each of the eight colors X, define qx{x) = dx{fx{x))- The level sets of qx are simply Cartesian 
products of the three dimensional level sets of dx with complementary nine-dimensional subspaces 
of M^^ , so they are convex and each qx is quasiconvex. 

It remains to formulate the requirement that our output gamut lie within the intersection of 
the input gamuts. If we are given n input gamuts, form a halfspace H^j (with < i < n and 
< J < 6) for each of the six facets of each of these parallelepipeds, and for each color X form 
a nested convex family = {x E M^^ | fx{x) G Hij} that ignores its argument A and 

returns a constant halfspace. We can then represent the problem of finding a feasible gamut that 
minimizes the maximum distance from one of its corners to the corner's ideal location as the 
quasiconvex program formed by the eight quasiconvex functions qx together with the 48n nested 
convex families Ki,j,x- 

4.6 Analysis of Backtracking Recurrences 

In this section we discuss another application of quasiconvex programming, in the automated anal- 
ysis of algorithms, from our paper [33]. There has been much research on exponential-time exact 
algorithms for problems that are NP-complete (so that no polynomial time solution is expected); 
see [6,12,22,29,30,32,44,75,82] for several recent papers in this area. Although other techniques 
are known, many of these algorithms use a form of backtracking search in which one repeatedly 
performs some case analysis to find an appropriate structure in the problem instance, and then uses 
that structure to split the problem into several smaller subproblems which are solved by recursive 
calls to the algorithm. 

For example, as part of a graph coloring algorithm [30] we used the following subroutine for 
listing all maximal independent sets of a graph G that have at most k vertices in the maximum 
independent set (we refer to such a set as a k-MiS). The subroutine consists of several different 
cases, and applies the first of the cases which is found to be present in the input graph G: 

— If G contains a vertex v of degree zero, recursively list each {k — 1)-MIS in G \ {-y} and append 
V to each listed set. 

— If G contains a vertex v of degree one, with neighbor ti, recursively list each {k — 1)-MIS in 
G \ N{u) and append u to each listed set. Then, recursively list each {k — 1)-MIS in G \ {u, v} 
and append v to each listed set. 

— If G contains a path vi-V2'V^ of degree-two vertices, then, first, recursively list each (k — 1)-MIS 
in G \ N{vi) and append vi to each listed set. Second, list each (k — 1)-MIS in G \ N{v2) and 
append V2 to each listed set. Finally, list each {k — 1)-MIS in G \ {{vi} U N^vs)) and append ^3 
to each listed set. Note that, in the last recursive call, vi may belong to N{vs) in which case 
the number of vertices is only reduced by three. 

— If G contains a vertex v of degree three or more, recursively list each fc-MIS in G \ {v}. Then, 
recursively list each [k — 1)-MIS in G \ N(v) and append v to each listed set. 

Clearly, at least one case is present in any nonempty graph, and it is not hard to verify that any 
/c-MIS will be generated by one of the recursive calls made from each case. Certain of the sets 
generated by this algorithm as described above may not be maximal, but if these non-maximal 
outputs cause difficulties they can be removed by an additional postprocessing step. We can bound 
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Table 1. A recurrence arising from unpublished work with J. Byskov on graph coloring algorithms, 
taken from [33]. 
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the worst-case number of output sets produced by this algorithm as the solution to the following 
recurrence in the variables n and k: 
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As base cases, T(0,0) = 1, T{n,—1) = 0, and T{n,k) = for k > n. Each term in the overall 
maximization of the recurrence comes from a case in the case analysis; the recurrence uses the 
maximum of these terms because, in a worst-case analysis, the algorithm has no control over which 
case will arise. Each summand in each term comes from a recursive subproblem called for that 
case. It turns out that, for the range of parameters of interest n/A < k < n/3, the recurrence above 
is dominated by its last two terms, and has the solution T{n,k) = (4/3)"' (3^/4^)'^. We can also 
find graphs having this many fc-MISs, so the analysis given by the recurrence is tight. Similar but 
somewhat more complicated multivariate recurrences have arisen in our algorithm for 3-coloring [29] 
with variables counting 3- and 4-value variables in a constraint satisfaction instance, and in our 
algorithm for the traveling salesman problem in cubic graphs [32] with variables counting vertices, 
unforced edges, forced edges, and 4-cycles of unforced edges. Another such recurrence, of greater 
complexity but with the same general form, is depicted in Tabled 

We would like to perform this type of analysis algorithmically: if we are given as input a 
recurrence such as the ones discussed above, can we efficiently determine its asymptotic solution, 
and determine which of the cases in the analysis are the critical ones for the performance of the 
backtracking algorithm that generated the recurrence? We showed [33] that these questions can be 
answered automatically by a quasiconvex programming algorithm, as follows. 

Let X denote a vector of arguments to the input recurrence, and for each term in the input 
recurrence define a univariate linear recurrence, by replacing x with a weighted linear combination 
^ = w ■ X throughout. For instance, in the /c-bounded maximal independent set recurrences, the 
four terms in the recurrence lead to four linear recurrences 



We can solve each of these linear recurrences to find constants q such that tj(^) = O(c^); it follows 
that, for any weight vector w, T{x) = 0(maxc^'^). 

This technique only yields a valid bound when each linear recurrence is solvable; that is, when 
each term on the right hand side of each linear recurrence has a strictly smaller argument than the 
term on the left hand side. In addition, different choices of w in this upper bound technique will 
give us different bounds. 

To get the tightest possible upper bound from this technique, for x = nt where t is a fixed target 
vector, constrain w -t = 1 (this is a normalizing condition since multiplying by a scalar does not 
affect the overall upper bound), and express q as a function Cj = qi{w) of the weight vector w; set 
Ci = +00 whenever the corresponding linear inequality has a right hand side term with argument 
greater than or equal to that on the left hand side. We show in [33] that these functions Qi are 
quasiconvex, as their level sets can be expressed by the formula 
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Fig. 14. The Tukey depth of the point marked with the + sign is three: there is a halfplane 
containing it and only three sample points (shown as solid disks); or, equivalently, three points can 
be removed from the sample set to place the test point outside the convex hull of the remaining 
points (shaded). 

where the right hand side describes a level set of a sum of convex functions of w. Therefore, we can 
find the vector w minimizing maxj gj (w) as a quasiconvex program. The value A of this quasiconvex 
program gives us an upper bound T{nt) = 0(A") on our input recurrence. 

In the same paper, we also show a lower bound T{nt) = f2{X""n~^), so the upper bound is tight to 
within a factor that is polylogarithmic compared to the overall solution. The lower bound technique 
involves relating the recurrence solution to the probability that a random walk in a certain infinite 
directed graph reaches the origin, where the sets of outgoing edges from each vertex in the graph 
are also determined randomly with probabilities determined from the gradients surrounding the 
optimal solution of the quasiconvex program for the upper bound. 

4.7 Robust Statistics 

If one has a set of n observations Xi S M, and wishes to summarize them by a single number, 
the average or mean is a common choice. However, it is sensitive to outliers: replacing a single 
observation by a value far from the mean can change the mean to an arbitrarily chosen value. 
In contrast, if one uses the median in place of the mean, at least n/2 observations need to be 
corrupted before the median can be changed to an arbitrary value; if fewer than n/2 observations 
are corrupted, the median will remain within the interval spanned by the uncorrupted values. In 
this sense, the median is robust while the mean is not. More generally, we define a statistic to be 
robust if its breakdown point (the number of observations that must be corrupted to cause it to 
take an arbitrary value) is at least cn for some constant c > 0. 

If one has observations Xi € M"^, it is again natural to attempt to summarize them by a single 
point X G M*^. In an attempt to generalize the success of the median in the one-dimensional prob- 
lem, statisticians have devised many notions of the depth of a point, from which we can define a 
generalized median as being the point of greatest depth [43,46,61-63,74,90,95]. Of these definitions, 
the most important and most commonly used is the Tukey depth [46,90], also known as half space 
depth or location depth. According to this definition, the depth of a point x (which need not be one 
of our sample points) is the minimum number of sample points contained in any halfspace that 
contains x (Figure E}- The Tukey median is any point of maximum depth. It follows by applying 
Helly's theorem to the system of halfspaces containing more than dn/{d + 1) observations that, 
for observations in W^, the Tukey median must have depth at least n/{d + 1). This depth is also 
its breakdown point, so the Tukey median is robust, and it has other useful statistical properties 
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as well, such as invariance under affine transformations and the ability to form a center-outward 
ordering of the observations based on their depths. 

There has been much research on the computation of Tukey medians, and of other points 
with high Tukey depth [14,19-21,51,57-59,65,73,79,87]. Improving on many previously published 
algorithms, Chan [14] found the best bound known for Tukey median construction, 0{n log n-\-n'^~^) 
randomized expected time, using his implicit quasiconvex programming technique. 

Let B he a bounding box of the sample points. Each d-tuple t of sample points that are in 
general position in R*^ defines a hyperplane that bounds two closed halfspaces, H^' and . If we 
associate with each such halfspace a number or that counts the number of sample points 
in the corresponding halfspace, then the pairs {B n H^, —kf) can be used to form a generalized 
longest intersecting prefix problem, as defined in Section 12.51 borrowing the terminology of LP- 
type problems, call any such pair a constraint. The solution to the quasiconvex program defined 
by this set of constraints is a pair (fe, x) where k is minimal and every halfspace with more than 
k samples contains x. If a halfspace H contains fewer than n — k samples, therefore, it does not 
contain x, so the depth of x is at least n — k. Any point of greater depth would lead to a better 
solution to the problem, so x must be a Tukey median of the samples, and we can express the 
problem of finding a Tukey median as a quasiconvex program. This program, however, has 0{n'^) 
constraints, larger than Chan's claimed time bound. To find Tukey medians more quickly, Chan 
applies his implicit quasiconvex programming technique: we need to be able to solve constant sized 
subproblems in constant time, solve decision problems efficiently, and partition large problems into 
smaller subproblems. 

It is tempting to perform the partition step as described after Theorem El by dividing the set 
of samples arbitrarily into d + 1 equal-sized subsets and using the complements of these subsets. 
However, this idea does not seem to work well for the Tukey median problem: the difficulty is that 
the numbers k^ do not depend only on the subset, but on the whole original set of sample points. 

Instead, Chan modifies the generalized longest intersecting prefix problem (in a way that doesn't 
change its optimal value) by including a constraint for every possible halfspace, not just those 
halfspaces bounded by d-tuples of samples. There are infinitely many such constraints but that 
will not be problematic as long as we can satisfy the requirements of the implicit quasiconvex 
programming technique. To perform the partition step for this technique, we use a standard tool 
for divide and conquer in geometric algorithms, known as e-cuttings. We form the projective dual of 
the sample points, which is an arrangement of hyperplanes in R"^; each possible constraint boundary 
is dual to a point in R'^ somewhere in this arrangement, and the number k^ for the constraint 
equals the number of arrangement hyperplanes above or below this dual point. We then partition 
the arrangement into a constant number of simplices, such that each simplex is crossed by at 
most en hyperplanes. For each simplex we form a subproblem, consisting of the sample points 
corresponding to hyperplanes that cross the simplex, together with a constant amount of extra 
information: the simplex itself and the numbers of hyperplanes that pass above and below it. Each 
such subproblem corresponds to a set of constraints dual to points in the simplex. When recursively 
dividing a subproblem already of this form into even smaller sub-subproblems, we intersect the 
sub-subproblem simplices with the subproblem simplex and partition the resulting polytopes into 
smaller simplices; this increases the number of sub-subproblems by a constant factor. In this way 
we fulfill the condition of Theorem El that we can divide a large problem into a constant number of 
subproblems, each described by an input of size a constant fraction of the original. 

Subproblems of constant size may be solved by constructing and searching the arrangement dual 
to the samples within the simplex defining the subproblem. It remains to describe how to perform 
the decision algorithm needed for Theorem O Decision algorithms for testing the Tukey depth of 
a point were already known [78,80], but here we need to solve a slightly more general problem due 
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to the extra information associated with each subproblem. Given k, x, and a subproblem of our 
overall problem, we must determine whether there exists a violated constraint; that is, a halfspace 
that is dual to a point in the simplex defined by the subproblem, and that contains more than k 
sample points but does not contain x. Let H be the hyperplane dual to x, and A be the simplex 
defining the subproblem. If there exists a violated constraint dual to a point h & A, we can assume 
without loss of generality that either h £ H oi h is on the boundary of A; for, if not, we could find 
another halfspace containing as many or more samples by moving h along a vertical line segment 
until it reaches either H or the boundary. Within H and each boundary plane of the simplex, we 
can construct the {d— l)-dimensional arrangement formed by intersecting this plane with the planes 
dual to the sample points, in time 0(nlogn + n'^~^). Within each face of these arrangements, all 
points are dual to halfspaces that contain the same number of samples, and as we move from face 
to face, the number of sample points contained in the halfspaces changes by ±1, so we can compute 
these numbers in constant time per face as we construct these arrangements. By searching all faces 
of these arrangements we can find a violated constraint, if one exists. 

To summarize, by applying the implicit quasiconvex programming technique of Theorem 
to a generalized longest intersecting prefix problem, using e-cuttings to partition problems into 
subproblems and {d — l)-dimensional arrangements to solve the decision algorithm as described 
above, Chan [14] shows how to find the Tukey median of any point set in randomized expected 
time 0(n log n + n'^~^). 

5 Conclusions 

We have introduced quasiconvex programming as a formalization for geometric optimization in- 
termediate in expressivity between linear and convex programming on the one hand, and LP-type 
problems on the other. Quasiconvex programs are capable of expressing a wide variety of geometric 
optimization problems and applications, but are still sufficiently concrete that they can be solved 
both by rapidly converging numeric local improvement techniques and (given the assumption of 
constant-time primitives for solving constant-sized subproblems) by strongly-polynomial combina- 
torial optimization algorithms. The power of this approach is demonstrated by the many and varied 
applications in which quasiconvex programming arises. 
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